
National Aeronautics and Space Administration 


Numerical Computation of a 
Continuous-Thrust State 
Transition Matrix Incorporating 
>\c jra'Q/H a rdware and 
Ephemeri Models 


Donald Ellison 
Bruce Conway 
Jacob Englander 


NAVIGATION & MISSION DESIGN BRANCH 

NASAGSFC 


www.nasa.gov 





■ Donald Ellison - Ph. D. Student, Department of Aerospace Engineering 
University of Illinois at Urbana-Champaign 

■ Bruce Conway - Professor, Department of Aerospace Engineering 
University of Illinois at Urbana-Champaign 

■ Jacob Englander - Aerospace Engineer, Navigation & Mission Design 
Branch, Code 595, NASA Goddard Space Flight Center 


NAVIGATION & MISSION DESIGN BRANCH, CODE 595 

NASA GSFC 


2 





■ Evolutionary Mission Trajectory Generator (EMTG) 

■ Finite-burn low-thrust (FBLT) transcription 

■ FBLT vs. multiple gravity-assist with low-thrust (MGALT) transcription 

■ FBLT match point constraint gradient calculation 

■ Continuous-thrust state transition matrix calculation 

■ Adaptive-step RK7(8)13M integrator 

■ FBLT match point time of flight gradient calculation 

■ Numerical Example 

■ Summary 
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The Evolutionar 
Generator (EMT 


■ Automated interplanetary low-thrust mission design tool 

■ Capabilities: 

- Global optimization of both the flyby sequence and the trajectory 

No initial guess required, global search capability provided by monotonic basin 
hopping (MBH) algorithm 

Works in multiple reference frames (interplanetary, moon tours, Earth orbit, etc) 

- Integration with SPICE ephemerides 

Up-to-date models of launch vehicles, thrusters, power systems 

- Easy to use graphical user interface 

■ The purpose of EMTG is to automate almost all of the work so that an analyst can 
investigate a wide range of mission options in very little human time 

■ Computer time is CHEAP. Analyst time is EXPENSIVE 
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J. A. Englander, D. H. Ellison, and B. A. Conway, “Global Optimization of Low-Thrust, Multiple-Flyby 
Trajectories at Medium and Medium-High Fidelity,” AAS/AIAA Space Flight Mechanics Meeting, Santa 
Fe, NM, January 26 - 30 2014. 
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FBLT vs. MG 


p 

Description 

Notes 

to 

Launch epoch 


Voo 

Magnitude of the outgoing velocity asymptote 

first phase only 

RA 

Right ascension of launch asymptote 

first phase only 

DEC 

Declination of launch asymptote 

first phase only 

TOF 

Phase time of flight 

N p 

rrif 

Phase final mass 

N p 

u 

Segment control vector 

one per FBLT segment (3 N x N p ) 

Isp 

Propulsion system specific impulse 

only for VSI-capable engines (N x N p ) 

Vooo 

Phase initial excess velocity vector 

all but the first phase (3N P — 1) 

Voo , 

Phase final excess velocity vector 

all phases, except the final phase of a rendezvous ( N p ) 


Table 1: Typical decision vector for a low-thrust mission FBLT mission. 
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FBLT as a No 


I 


n: Constraints 


MGALT and FBLT share identical decision and constraint vectors 

Control bounds, flyby model and TOF constraint gradients may be 
provided fully analytically 
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i: Constraints 


■ FBLT transcription results in a large, sparse NLP 

■ We want to calculate the gradient of each NLP constraint w.r.t. each NLP parameter 
(KKT first order necessary conditions) — ► form the problem Jacobian 

■ SNOPT will do this with finite differencing... this is EXPENSIVE 

■ FBLT is already much slower than MGALT due to numerical integration vs. analytic 
Kepler, so it is even more important to supply user-defined derivatives 

■ The most challenging are the match point continuity constraints: state transition 
matrices 




Op Op Op 
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Match Point 
Calculation 


■ Traditional six-by-six STM is augmented to include mass and the 
current time step’s control parameters u xt u y and u z 

■ Match point constraint vector sensitivities w.r.t. controls make up the 
majority of the dense entries in the Jacobian 

■ An STM chain is constructed beginning with the segment of interest and 
proceeding using “stripped” STMs to the match point 
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Calculation 


■ Unlike the two-body problem (see Battin for example), analytic expressions for the 
perturbation STM entries don’t exist for true continuous thrust 

■ We can solve for them numerically 

■ Differential equations for the STM entries are appended to the physics EOM’s and 
integrated alongside them 

■ We must calculate the state propagation matrix (A) at each integration time step 

■ A also referred to as the fundamental solution/set 

3> — A3> 

= f X y z v x Vy v z m 1 T A — dX. 

A ~ dX 

* r . # i t 

— f — | i* r 7TI %l x Uy U z S 12 * * * Sioio 

S. P. Hughes, D. S. Cooley, and J. J. Guzman, “A Direct Method for Fuel Optimal Maneuvers of Dis- 
tributed Spacecraft in Multiple Flight Regimes,” Proceedings of the AAS/AIAA Space Flight Mechanics 
Meeting , Copper Mountain, Colorado , No. AAS-05-158, January 23-27 2005. 

S. R Hughes and E. Dove, “Optimal Control and Near Optimal Guidance for the Magnetospheric Multi- 
Scale Mission (MMS),” Proceedings of the AAS/AIAA Astrodynamics Specialist Conference , Pittsburgh , 

Pennsylvania , No. AAS-09-330, August 9-13 2009. 



FBLT State Tr 
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Calculation 
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Calculation 
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Calculation 
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Calculation 
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Launch Vehi 
Modeling 



■ Launch vehicles are modeled using a polynomial fit 

m delivered = (l — + e LV^ 3 + d LV C f + C LV C 2 + b LV C3 + 

where o LV is launch vehicle margin and C 3 is hyperbolic excess velocity 

■ Thrusters are modeled using either a polynomial fit to published thrust and mass flow rate 
data 

ril m ax = &F P 4 T dp P 2 + Cp P 2 + bp P + Cip 
Tmax = e T P 4 + d T P 3 + CpP 2 + b T P + dp 


or, when detailed performance data is unavailable 


max 
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max 


max 
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Power is modeled by a standard polynomial model 

v 4-11 + 12 , 

Yo + r + r 2 


1 + y 3 r + y 4 r 2 


(1 - t ) 1 


P 0 is the power at beginning of life at 1 AU, r is the solar array degradation constant and t is 
the time since launch in years 
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EMTG’s Integra 


■ 


>)13M 


Two main components 


- Adaptive step sizing algorithm 

- 8 th order embedded explicit Runge-Kutta step calculation method 

■ Local truncation error between 7 th and 8 th order solutions is used to adaptively 
size the integration step 

^0^1 . . . tn ^n+l ... tf 

■ Each FBLT 
segment h N is 
broken into 8 t 
sub-steps 

Hn = ^ Si 


II I 

do Si ■ • • s„ 

hx 
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■ After each sub-step, the local truncation error between the 7 th and 8 th order 
solutions is calculated (Ae) 

■ If Ae > AT then S = 0.98, the time step will be reduced and the sub-step is 
reintegrated 

■ else S = 1.01, and the time step length is increased, potentially saving 
computation time 

■ q is the order of the relative error of the method, q = 8 for RK7(8)13M 

■ The integration proceeds until integration across the full FBLT time step has 
been completed to within the user-specified error tolerance 



A. M. Ghosh, Multi-Cubesat Mission Planning Enabled Through Parallel Computing. PhD thesis, 
University of Illinois at Urbana-Champaign, May 2013. 
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EMTG’s Integra 


■ 


7)1 3M 


Calculate the gradient and step at the left hand side of the RK sub-step 

ki = S n ■ f (t n , X, 

At each of the 13 stages, calculate the state sample point and the 
gradient/step at that point... store the aradient information 
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4. Update time and adjust the step size ^ ^ ^ ^ 
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P. J. Prince and J. R. Dormand, “High order embedded Runge-Kulta formulae,” Journal of Computa- 
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P. J. Prince and J. R. Dormand, “High order embedded Runge-Kulta formulae,” Journal of Computa- 
tional and Applied Mathematics, Vol. 7, No. 1, 1981, pp. 67-75. 
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Match Point Ti 


ivatives 


■ If we had an analytic equations for propagating the trajectory, TOF derivatives 
would be fairly simple (to first order) 

5r rn = v m * STOF 


■ Unfortunately, we have a numerical approximation, so computing TOF 
derivatives involves taking derivatives of the components of the integrator itself 

■ Derivatives of both the adaptive step routine as well as the RK step calculation 
must be accumulated as the integration is performed 

0X n+ i = 8X n ™ dki 
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* ■ 
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dtcb-B/c dT 


TOF Derivati 


■ State gradient TOF derivatives for the spacecraft in turn require sensitivity of 3 rd 
body position to changes in TOF 

■ EMTG uses the SPICE kernels to calculate solar system body states 

■ Analytical derivatives are only possible if you can access SPICE internally 

■ SPICE uses Chebyshev polynomial interpolation for planetary positions 

■ Fortunately, for planets only, the velocity polynomial is the time derivative of the 
position polynomial (we have direct access to the position curve derivative) 

■ This is not true for planetary satellites, velocity curve is fit separately 

Another ephemeris reader (FIRE, Arora and Russell, 2008) would solve this 

■ EMTG’s power solar electric hardware models are TOF dependent 
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TOF Derivati 


■ State gradient TOF derivatives for the spacecraft in turn require sensitivity of 3 rd 
body position to changes in TOF 

■ EMTG uses the SPICE kernels to calculate solar system body states 

■ Analytical derivatives are only possible if you can access SPICE internally 

■ SPICE uses Chebyshev polynomial interpolation for planetary positions 

■ Fortunately, for planets only, the velocity polynomial is the time derivative of the 
position polynomial (we have direct access to the position curve derivative) 

■ This is not true for planetary satellites, velocity curve is fit separately 

Another ephemeris reader (FIRE, Arora and Russell, 2008) would solve this 

■ EMTG’s power solar electric hardware models are TOF dependent 
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^■3x3 

d r c b-3b 

3b— s/c 

r 3 

# 3b— s/c _ 
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TOF Derivati 


■ State gradient TOF derivatives for the spacecraft in turn require sensitivity of 3 rd 
body position to changes in TOF 

■ EMTG uses the SPICE kernels to calculate solar system body states 

■ Analytical derivatives are only possible if you can access SPICE internally 

■ SPICE uses Chebyshev polynomial interpolation for planetary positions 

■ Fortunately, for planets only, the velocity polynomial is the time derivative of the 
position polynomial (we have direct access to the position curve derivative) 

■ This is not true for planetary satellites, velocity curve is fit separately 

Another ephemeris reader (FIRE, Arora and Russell, 2008) would solve this 

■ EMTG’s power solar electric hardware models are TOF dependent 


®^cb— sjc 

dTOF 


dv cb—s/c 
d^cb— sjc 


d^cb—s/c ^^ch—s/c 

dTOF dr cb . 3b 


pr cb - 


3 b 


dTOF 


+ 


dVrh- 


ch—sj C 


dm 


sf C 


dm s/c s/c 

DTOF + 9T 


dT 

OTOF 
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TOF Derivati 


■ State gradient TOF derivatives for the spacecraft in turn require sensitivity of 3 rd 
body position to changes in TOF 

■ EMTG uses the SPICE kernels to calculate solar system body states 

■ Analytical derivatives are only possible if you can access SPICE internally 

■ SPICE uses Chebyshev polynomial interpolation for planetary positions 

■ Fortunately, for planets only, the velocity polynomial is the time derivative of the 
position polynomial (we have direct access to the position curve derivative) 

■ This is not true for planetary satellites, velocity curve is fit separately 

Another ephemeris reader (FIRE, Arora and Russell, 2008) would solve this 

■ EMTG’s power solar electric hardware models are TOF dependent 


®^cb— sjc 

dTOF 


di* cb—s/c dr cb—s/c s/c 

dr cb _ s / c dTOF dr c& _ 36 


d^ch—3b cb—s/c 

DTOF + dm g/c 


dui s / c s / c 

OTOF + 9T 


dT 

1 WOF 
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Example: Sol 

M ! ic Earth to Jupiter 


Parameter 

Value 

Initial mass 

3000 kg 

Earliest allowed launch date 

January 1 st, 2015 

Latest allowed launch date 

unbounded 

Initial Voo 

up to 6.97 km/s 

Maximum flight time 

unbounded 

Arrival type 

low-thrust rendezvous 

Thruster 

NEXT - high thrust variant 

Thrust coefficients 

gt = 0.09591 d T = —1.98537 c T - 11.47980 b T = 15.0G977 a T - 14.51552 

Thruster duty cycle 

0.9 

Mass flow rate coefficients 

e F = 0.01492 d F = -0.27539 c F = 1.60966 b F = -2.53056 a F = 3.22089 

Solar power coefficients 

70 = 1.32077 71 = -0.10848 72 = -0.11665 73 = 0.10843 74 = -0.01279 

Solver parameters 


Flyby sequence 

Earth-Jupiter direct (E-J) 

Number of time -steps 

40 

Ephemeris 

SPICE 

SNOPT feasibility tolerance 

1.0c-5 

Objective function 

maximize final mass 

Number of NLP parameters 

126 

N li mbe r of c on strai n ts (i nc 1 udi n g obj ect i ve fu n ct ion) 4 8 

Dense Jacobian entries 

1002 

Table 2: Earth to Jupiter low-thrust direct transfer problem assumptions 
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Example: Sol 

trie Earth to Jupiter 


x (km) 


5 i 



5 



-5 




NAVIGATION & MISSION DESIGN BRANCH, CODE 595 

NASAGSFC 


27 


z (km) le8 




Example: Sol 

trie Earth to Jupiter 



NAVIGATION & MISSION DESIGN BRANCH, CODE 595 

NASAGSFC 


28 




Example: Sol 

trie Earth to Jupiter 
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Example: Sol 



Jupiter 


■ SNOPT was allowed to execute for 100 major iterations 

■ STM-based Jacobian run converged after 80 majors 

■ Finite differenced Jacobian run hit the 100 major iteration limit 

■ Both are equally effective at solving the problem from an MGALT initial 
guess 

■ Numerically computed STMs afford 13-15 times execution speed increase 


Metric 

STM computed Jacobian 

Finite Differenced Jacobian 

Final mass delivered to Jupiter 

2117.68 kg 

21 17.53 kg 

SNOPT execution time 

44.04 s 

672.91 s 


Table 3: STM vs. finite differenced Jacobian SNOPT performance metrics 
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Next Steps 


■ Low-thrust preliminary design cadence has been drastically increased 

■ EMTG-General Mission Analysis Tool (GMAT) work flow has been 
improved 

■ Although GMAT can converge from an MGALT or FBLT initial guess, 
for complex problems sometimes even an FBLT cannot converge from 
MGALT 


For cases when MGALT is insufficient, it is crucial that initial guesses 
for GMAT be generated as fast as possible using FBLT 

Launch vehicle, power and thruster models are currently being 
integrated with GMAT 

It is recommended that the STM described in this presentation be 
incorporated into GMAT’s solvers as it can handle arbitrarily complex 
dynamics (i.e. SRP, gravitational harmonics could be added relatively 
easily) 

EMTG and GMAT native integrators are also quite similar 



NAVIGATION & MISSION DESIGN BRANCH, CODE 595 

NASAGSFC 


Conclusions 


■ The state transition matrix increases the efficiency of the finite-burn low-thrust 
transcription and allows for the tempo at which low-thrust preliminary design is 
performed to be increased, saving analysts’ time. 

■ It importantly includes accurate hardware models that can have a major influence 
on mission feasibility even at the preliminary design stage 

■ This work enables FBLT to be used in the framework of a global optimizer for 
many problems 

■ With EMTG’s performance increases serving as an example, the logical next step 
would be to repeat the same process inside the GMAT high-fidelity tool 

■ This work has extensions in other areas of interest: 

- Spacecraft guidance (method of adjoints and the guidance matrix) 
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Thank you 


This work is supported by Rich Burns and was completed 
under contract with NASA GSFC NNX14AD27G 

EMTG is available open-source at: 
https ://so u reef o rg e . n et/p ro j ects/e m to/ 
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